Stretched exponentials in the dynamics and at the ground state of the plaquette Ising 

model 
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Simulations show that the autocorrelation function C(t) in the d — 3 Ising model with a pla- 
quette interaction has a stretched-exponential decay in a supercooled liquid phase. Such a decay 
characterizes also some ground-state probability distributions obtained from the numerically exact 
counting of up to 10 450 configurations. A related model with a strongly degenerate ground state 
but lacking glassy features does not exhibit such a decay. Tensionless modification of the droplet 
model might explain stretched-exponential decay of C(t) even in three-dimensional systems. 



Assuming that liquids constitute a homogeneous and 
continuum medium, one can explain many of its prop- 
erties such as viscosity, diffusion or chemical reaction 
rates. However, cooling below melting point causes dras- 
tic changes in the dynamics and a homogeneous approach 
is no longer legitimate One of the landmarks of 

this supercooled regime is a slower than exponential de- 
cay of autocorrelation functions. Although this decay 
is often fitted with the so-called stretched exponentials 
exp[— {t/r)P], it is sometimes questioned as being merely 
a phenomenological fit to experimental data without any 
clear microscopic mechanism Q. 

It is becoming commonly accepted that supercooled 
liquids are dominated by dynamically generated hetero- 
geneities. They have different sizes and lifetimes and 
their dynamics is thus very complex. It would be de- 
sirable to relate this dynamics with some kind of an 
Ising model the dynamics of which is relatively under- 
stood. Particularly interesting in this context might 
be the droplet model, which predicts in some cases the 
stretched-exponential decay of autocorrelation functions 
According to the droplet model, however, such a de- 
cay should hold only for low-dimensional Ising systems, 
and in the most interesting three-dimensional case, an 
exponential decay is expected. 

A promising statistical mechanics approach to glassy 
systems refers to lattice models. Various models 
of glasses, including kinetically constrained or spin- 
facilitated ones have already been examined Q. Al- 
though they exhibit an interesting slow dynamics, their 
thermodynamics is very often trivial. A more comprehen- 
sive description of glasses might be thus sought among 
Ising models. 

In the present Letter, we examine the Ising model with 
plaquette interaction, that is defined using the following 
Hamiltonian 



H 



Si Sj Sk Si , 



±1, 



1,2,. ..,L d (1) 



where summation is over elementary plaquettes of the d- 
dimensional Cartesian lattice of the linear size L with pe- 
riodic boundary conditions. For d — 3, model ((T|) shares 
a number of properties with glassy systems. In partic- 



ular, it exhibits a strong metastability [H, 0] and a very 
slow (perhaps logarithmically slow) coarsening dynamics 
7}. Moreover, aging 

6] 

and cooling-rate effects [8| are 
consistent with expectations for ordinary glasses. Let us 
also notice that even the d = 2 version of the model, 
despite trivial thermodynamics, exhibits an interesting 
glassy behaviour 0, ■ 

To examine the dynamics of model (TTJ) in the super- 
cooled liquid phase, we calculated the spin-spin autocor- 
relation function that is defined as 



C(t) 



(2) 



Using a standard Metropolis dynamics, we simulated 
the model for the d = 3 case at several temperatures 
and measured C(t). The results of these simulations 
are shown in FigjTJ Fitting our data to the function 
aexp [— (i/r)' 3 ], we obtained the exponent (3, and its tem- 
perature dependence is shown in the inset. Let us notice 
that earlier studies of model ([T]) locate the glassy tran- 
sition close to T = T g s» 3.4 and the first-order melt- 
ing transition (located on the comparison of free energies 
and simulations with nonhomogeneous initial conditions) 
around T — T m w 3.6 0, @, Hlj . One can notice that for 
T < T m an appreciable departure from the exponential 
(/3 = 1) decay is seen. It is not clear to us whether (3 
becomes smaller than 1 precisely at T m or at a some- 
what larger value, as our data might suggest. Studies of 
long-time evolution of glassy systems using molecular dy- 
namics simulations of realistic systems are computation- 
ally very demanding, but in a model with a controlled 
frustration, Shintani reported a similar temperature de- 
pendence of /3 12;]. A stretched-exponential behaviour 
of the energy autocorrelation function with j3 decreasing 
upon approaching T g has already been reported for model 
([T]) However, energy is a global variable and it is not 
clear whether such a quantity correctly probes the het- 
erogenoeus dynamics of supercooled liquids. There are 
also some reports of a stretched exponential behaviour 
in other Ising-like lattice models with glassy features [l3[ 
but it might be a consequence of a reduced (d = 2) di- 
mensionality, which according to the droplet model Q 
might imply such a decay of C(t). 



FIG. 1. The time dependence of the spin-spin autocorre- 
lation function C(t) for (from above) T = 3.43, 3.45, 3.5 
, 3.6, 3.8, 4.0, and 4.2. At each temperature, simulations 
start from a random initial configuration, which is relaxed 
for a certain time (~ 10 3 ), and then measurements are made 
during t ~ 10 7 (L — 50). A unit of time corresponds to a 
single on average update per spin. The inset shows the ex- 
ponent P obained from the fit to the stretched exponential 
(aexp — (t/r)' 9 ). The arrow indicates location of the first- 
order melting transition in the model. 



As we show below, stretched exponentials appear in 
model (QJ also in a much different context. Namely, they 
characterize some ground-state probability distributions. 
A calculation of these distributions is possible because a 
strongly degenerate ground state has in fact a relatively 
simple structure. For example, for d = 2 all its config- 
urations can be obtained from a reference configuration 
(e.g., all spins +) by flipping vertical and horizontal rows 
of spins (FigHJ), which leads to its 2 2i_1 dc generacy 0. 

To illustrate the calculations, let us examine the 
susceptibility-like variable xi = T^^2ij^i^j = 
jY(J2i Si) 2 ■ One can notice that to calculate X2 pre- 
cise distribution of these flipped rows is not needed, and 
it is sufficient to know only their numbers k\ and k 2 . In- 
deed, flipping fci horizontal rows, we reduce the number 
of + spins to N + = L 2 — k\L. The subsequent flip of 
k 2 vertical rows leads to N + = L 2 — (k\ + k2)L + 2k\k 2 . 
Using = — L 2 , we obtain that for (k\,k 2 ) 

configurations \2 — j?[L 2 — 2(fci + k 2 )L + 4kik 2 ] 2 . Of 
course, calculating the probability distribution of X2, one 
should take into account the multiplicity factor equal to 

( fei ) ( ki )■ 

The above considerations can be easily generalized to 
the d = 3 version. In this case, ground-state configura- 
tions can be obtained from the reference configuration 
by flipping entire two-dimensional planes. To character- 
ize a given ground-state configuration, we now need three 



FIG. 2. An example of a ground-state configuration in the 
d — 2 case obtained from the ferromagnetic state (all +) by 
flipping fci = 4 horizontal and &2 = 3 vertical rows of spins. 



numbers (fci, k 2 , k 3 ), for which we obtain 
X3 = -^{L 3 -2(k 1 + k 2 + k 3 )L 2 

+4(fcifc 2 + fcifc 3 + k 2 k 3 )L - 8fcifc 2 A;3] 



(3) 



The corresponding multiplicity factor equals 

( fci ) ( hi ) ( k 3 ) ' 

For a further analysis of the probability distribution, 
we resort to numerical caclulations. For d = 3 and a given 
system size L, we generate all (L— l) 3 triples (k\,k 2 , k 3 ), 
and using Eq.Q, we calculate Xs an d the corresponding 
multiplicity factor. Collecting the data in some bins, we 
obtain the required probability distribution P(% 3 ). Let 
us notice that computational complexity of generation of 
such triples is rather modest (~ L 3 ), which allows us to 
examine large systems (L ~ 500) within few seconds of 
CPU time. 

The results presented in Figj3]show that P(x3) has a 
maximum at x = 0. Let us notice, however, that the av- 
erage over all ground-state configurations (X3) = 1. This 
is a consequence of the flipping symmetry of the Hamil- 
tonian |T]), which implies that for any i ^ j the corre- 
sponding correlation function (SiSj) vanishes @. The 
only nonvanishing contribution to susceptibility comes 
from the case i = j and that implies (x 3 ) = 1. 

The above symmetry arguments are valid also at finite 
(and sufficiently high) temperature and recent Johnston's 
Monte Carlo simulations report indeed (X3) ~ 1 1A \- 



These simulations also show that at low temperature the 
susceptibility drops almost to zero. Our results (FigEJ) 
shed some light on such a finding: Monte Carlo simu- 
lations at low temperature select randomly one of the 
ground states (toward which the system slowly evolves) 
and since X3 — is the most probable value in P(x3), 
this is the value that is typically measured in Monte Carlo 
simulations. Let us notice that due to the strong degener- 
acy of the ground state, it is difficult to find the order pa- 
rameter that would distinguish low and high temperature 
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FIG. 3. The probability distribution for the susceptibility 
X3- The fact that X3 — is the most typical value of these 
distributions agrees with the recent Monte Carlo simulations 
for lattices of similar size showing that at low temperature the 
susceptibility drops to zero [lj. The inset shows P(x3 = 0) 
as a function of 1/L for L = 6,10, ... , 500. 
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FIG. 4. Logarithm of the probability distribution P(x3) f° r 
the three-dimensional version. For L — 500, the data nearly 
overlap with the least-square fit with j3 = 1/3. The inset 
confirms the stretched exponential decay with ft = 1/3 for 
nearly 120 decades. 



phases of the model {l}. Johnston's result suggests that 
the susceptibility might serve as such, but P(x = 0) for 
increasing L does not converge to unity (inset in FigJ3]) 
and there is a finite (albeit small) probability that simu- 
lations will select the ground state with x > 0. 

Further analysis of P(x3) (FigHJ) shows a slower than 

1/3 

exponential decay for large X3- Plotting against x 3 
shows an excellent linearity of the data for nearly 120 
decades and indicates that asymptotically P(X3) decays 

1/3 

as stretched exponential cxp(— a% 3 ). Let us also notice 
that for L = 500 the obtained probability distribution is 
based on the (numerically) exact counting of 2 3 ' 500 ~ 2 ~ 
10 450 ground-state configurations. 
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FIG. 5. Logarithm of the probability distribution P(xa) for 
the two-dimensional version. For L — 2000, the data nearly 
overlap with the least-square fit with /3 = 1/2. The inset 
confirms the stretched exponential decay with ft = 1/2 for 
nearly 120 decades. 



Although the d = 2 version of the model |T]) has a triv- 
ial thermodynamic behaviour, it exhibits an interesting 
glassy- like dynamic behaviour [§]. Calculating the dis- 
tribution P(x2), we also find the stretched exponential 
decay but with the exponent (3 = 1/2 (FigJSJ). 

It is not clear to us whether stretched exponential 
(static) distributions that we found at the ground state 
of model |T]) are related at all with its finite-temperature 
dynamic behaviour. However, a curious example comes 
from a more general version of model |T|) known as the 
gonihedric model. This model, introduced in the context 
of the discretized string theory [l5|], is described by the 
following Hamiltonian 



H — — 2k ^ ^ SiSj 



2 X! St S i 



1 — K 



SiSjSkSi 

(i,j,k,l) 

(4) 

where the first and the second summations are over the 
nearest and the next-nearest neighbours, respectively. 
For k = 0, the model reduces to the plaquette model 
([T]). A quite different thermodynamic and dynamic be- 
haviour is reported for n @, EH- In such a case, 
there is no metastability upon temperature changes and 
the dynamics does not exhibit glassy features. Also the 
flipping symmetry of the model @ is lower than in the 
case of k = 0. In particular, flipped planes or rows of 
spins cannot cross [16]. This simplifies the calculations 
we made for the plaquette model since now only one of 
the numbers ki, or k% might be nonzero. In the d = 3 
version of the gonihedric model for a configuration with 
k planes flipped, we thus obtain 



X3« 



-[L 3 - 2kL 2 } 2 = L[L - 2k} 2 . 



(5) 



and the multiplicity factor being equal to 3(^) (with 
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FIG. 6. Rescaled logarithm of the probability distribution 
P(X3k) for the three-dimensional version of the gonihedric 
model for k ^ 0. For L — 1000, the data nearly overlap 
with the straight line corresponding to exp (—x) confirming 
thus the asymptotic estimation JJJ. 



prefactor 3 corresponding to the number of directions of 
flipped planes). Using the above equation, one obtains 



() 



3-2^ \±(L- VP) 



(6) 



where 3 • 2 L is the degeneracy of the ground state. Using 

)(£/2)l(L/2)l 



the identity 

form (l/ 2 ) ? 

that for 1 < \/Wi < i 



(L-fc)ifc! — an< ^ the as y m pt°ti c 

after some calculations one obtains 



2exp(- 



2L 2 



'ttL 



(7) 



Our numerical data are in a very good agreement 
with the above estimation (FigJB]). Let us notice that 
P(X3k) has an exponential decay that in the thermody- 
namic limit L — > oo flattens and P(x3k) — ► 0. Thus 
the ground-state distribution in the gonihedric model for 
k =/= has a much different form than in the k — case. 
Perhaps it is only a coincidence that in the latter case, 
where the model exhibits glassy behaviour, the distribu- 
tion P(X3) nas a stretched-exponential behaviour. But 
one cannot exclude a more profound relation between 
zero-temperature statics and finite-temperature dynam- 
ics of model JT]). 

In the final part of our Letter, we would like to return 
to the problem of the decay of C(t) in the supercooled 
regime. It would be desirable to explain such a dynam- 
ics in terms of some kind of the droplet model 3]. The 
droplet model most likely provides a qualitatively correct 
description of the Ising dynamics 17] but it predicts the 
exponential decay of C(t) for d — 3 systems. It seems to 
us, however, that trying to use the droplet model to ex- 
plain the dynamics of model (fTJ , we might have to modify 
some of its assumptions. Indeed, the excess energy ei of a 



is proportional to its surfa ce ( e; ~ / 



d-l\ 



)del © 



might scale as ei ~ l d ~ 2 (jl floj. Actually, droplets with 
energy proportional to are also possible in model 
([T]) but in our opinion the long-term dynamics might be 
under the influence of mainly the (low-energy) tension- 
less droplets. Assuming e; ~ I 1 and repeating Lifshitz 
reasoning, one obtains that the deterministic motion of 
a droplet satisfies dl/dt ~ l~ 2 and the resulting lifetime 
of a domain of the initial size I thus scales as I 3 . Con- 
sequently, the form of the Boltzmann factor leads to the 
estimation C(t) ~ exp [— (i/r)^ -2 '/ 3 ]. Thus, tensionless 
domains in model ([1]) might imply stretched exponential 
behaviour in three-dimensional systems. 

The dynamics of model ([1]) in the supercooled liquid 
phase is very complex and most likely the entire spectrum 
of droplets is present. Droplets with energies scaling as 
l d ~ 2 and Z d 1 are only limiting cases and perhaps a more 
adequate description could be obtained asssuming e; ~ 
l^, where d— 2 < <f) < d— 1. In such a case, we are lead to 
C(t) - exp[-(i/T) 0/(d+1 ~«]. In our opinion, it would 
be desirable to calculate surface tension of droplets in 
model |T]) possibly confirming their tensionless nature. 

I thank Des Johnston for sending me his paper prior 
to the publication. 
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